library(tidyverse)
library(here)
library(phyloseq)
library(vegan)
theme_set(theme_bw())
filt_fun <- function (x, min_reads = 100, min_samples = 5) {
(sum(x) > min_reads) & (sum(x > 0) > min_samples)
}
ps <- readRDS(here('data','following_study','ps_rarefied.rds')) %>%
filter_taxa(function (x) sum(x) > 0, prune = TRUE)
sample_data(ps)[is.na(sample_data(ps))] <- 'NA'
# transform data into proportion
ps.prop <- ps %>% transform_sample_counts(function(x) x/sum(x))
# double check proportion transform
sample_sums(ps.prop) %>% summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
sam <- data.frame(sample_data(ps))
max.core <- parallel::detectCores()
alpha.diversity <- data.frame(sample_data(ps),estimate_richness(ps))
ggplot(alpha.diversity,aes(x = as.numeric(Household), y = Shannon, group = Household)) +
geom_point() + geom_line() + xlab('Household')
anova(lm(Shannon~Household, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.18604 1.5916 0.05425 .
## Residuals 49 5.7277 0.11689
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity,aes(x = Epileptic.or.Control, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
Paired t test
epileptic <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
control <- alpha.diversity %>% filter(Epileptic.or.Control == 'Control')
t.test(epileptic$Shannon, control$Shannon, paired = TRUE)
##
## Paired t-test
##
## data: epileptic$Shannon and control$Shannon
## t = -0.22089, df = 48, p-value = 0.8261
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.1556564 0.1248409
## sample estimates:
## mean difference
## -0.01540773
ANOVA with Household effect added
anova(lm(Shannon~Household + Epileptic.or.Control, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.186045 1.5607 0.06329 .
## Epileptic.or.Control 1 0.0058 0.005816 0.0488 0.82612
## Residuals 48 5.7219 0.119206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity) +
geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
facet_wrap(~Epileptic.or.Control) +
theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())
# here breed group with na is removed
anova(lm(Shannon~Household + Breed.Group..1., data = alpha.diversity %>% filter(Breed.Group..1. != 'NA')))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 45 7.4046 0.16455 1.5646 0.07819 .
## Breed.Group..1. 4 1.2267 0.30668 2.9160 0.03339 *
## Residuals 39 4.1016 0.10517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alpha.diversity.epi <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(alpha.diversity.epi, aes(x = Pheno.Y.N, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
t.test(Shannon~Pheno.Y.N, data = alpha.diversity.epi)
##
## Welch Two Sample t-test
##
## data: Shannon by Pheno.Y.N
## t = 0.98407, df = 33.649, p-value = 0.3321
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
## -0.1316909 0.3787797
## sample estimates:
## mean in group No mean in group Yes
## 2.856864 2.733320
ggplot(alpha.diversity, aes(x = Sex, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova(lm(Shannon~Household + Sex, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.18604 1.6309 0.04673 *
## Sex 1 0.2522 0.25218 2.2107 0.14360
## Residuals 48 5.4755 0.11407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(Shannon~Sex, data = alpha.diversity)
##
## Welch Two Sample t-test
##
## data: Shannon by Sex
## t = -0.83555, df = 88.411, p-value = 0.4057
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -0.2223635 0.0907199
## sample estimates:
## mean in group F mean in group M
## 2.792318 2.858140
ggplot(alpha.diversity, aes(x = as.numeric(Age..months.)/12, y = Shannon)) +
geom_point() + xlab('Age')
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
lm(Shannon~as.numeric(Age..months.), data = alpha.diversity) %>% summary()
## Warning in eval(predvars, data, env): NAs introduced by coercion
##
## Call:
## lm(formula = Shannon ~ as.numeric(Age..months.), data = alpha.diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90685 -0.30999 0.01342 0.28351 0.75221
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.745594 0.091949 29.860 <2e-16 ***
## as.numeric(Age..months.) 0.000838 0.001114 0.752 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3838 on 94 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.005986, Adjusted R-squared: -0.004588
## F-statistic: 0.5661 on 1 and 94 DF, p-value: 0.4537
plot_ord <- function(data, factor, method, distance){
data.ord <- ordinate(data, method = method, distance = distance)
p <- plot_ordination(data, data.ord, color = factor)
p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
model <- as.formula(paste0('dist.matrix~', formula))
if (!is_null(strata)) {strata <- df[,strata]}
result <- adonis2(model,
data = df,
permutations=permutations,
strata = strata,
parallel = core)
return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
beta.disp <- betadisper(dist.matrix, group = df[,group])
result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
return(result)
}
permanova(ps.prop, 'Household', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.2596 9.999e-05 ***
## Residual 49 5.5192 0.31119
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_ord(ps.prop, 'Epileptic.or.Control','MDS','bray')
plot_ord(ps.prop, 'Epileptic.or.Control','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2365 -- in the phylogenetic tree in the data you provided.
plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','bray')
## Run 0 stress 0.2086595
## Run 1 stress 0.2039329
## ... New best solution
## ... Procrustes: rmse 0.0483378 max resid 0.2803026
## Run 2 stress 0.2047374
## Run 3 stress 0.2280233
## Run 4 stress 0.2111629
## Run 5 stress 0.2389114
## Run 6 stress 0.2045404
## Run 7 stress 0.2093663
## Run 8 stress 0.2031223
## ... New best solution
## ... Procrustes: rmse 0.01902768 max resid 0.1332267
## Run 9 stress 0.224303
## Run 10 stress 0.2237001
## Run 11 stress 0.2042402
## Run 12 stress 0.222362
## Run 13 stress 0.2216106
## Run 14 stress 0.2174793
## Run 15 stress 0.2030407
## ... New best solution
## ... Procrustes: rmse 0.01326296 max resid 0.08083261
## Run 16 stress 0.2047467
## Run 17 stress 0.2053419
## Run 18 stress 0.2041372
## Run 19 stress 0.2120025
## Run 20 stress 0.2029529
## ... New best solution
## ... Procrustes: rmse 0.004458111 max resid 0.03007541
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 5: no. of iterations >= maxit
## 15: stress ratio > sratmax
plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2253 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1764419
## Run 1 stress 0.186396
## Run 2 stress 0.1762436
## ... New best solution
## ... Procrustes: rmse 0.02963319 max resid 0.1750311
## Run 3 stress 0.1853109
## Run 4 stress 0.1938895
## Run 5 stress 0.1849004
## Run 6 stress 0.1849142
## Run 7 stress 0.1821115
## Run 8 stress 0.1836205
## Run 9 stress 0.1795957
## Run 10 stress 0.1920684
## Run 11 stress 0.1777221
## Run 12 stress 0.1773028
## Run 13 stress 0.1834667
## Run 14 stress 0.1821115
## Run 15 stress 0.1861628
## Run 16 stress 0.1851304
## Run 17 stress 0.1814618
## Run 18 stress 0.1928126
## Run 19 stress 0.187659
## Run 20 stress 0.1743156
## ... New best solution
## ... Procrustes: rmse 0.02428976 max resid 0.1780691
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmax
permanova(ps.prop, 'Epileptic.or.Control', 'bray', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.1554 0.00876 0.8487 0.1022
## Residual 96 17.5803 0.99124
## Total 97 17.7358 1.00000
permanova(ps.prop, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV351 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.0546 0.01444 1.4069 0.0424 *
## Residual 96 3.7283 0.98556
## Total 97 3.7830 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.prop, 'Epileptic.or.Control', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00167 0.0016733 0.1117 10000 0.7339
## Residuals 96 1.43851 0.0149844
permdisp(ps.prop, 'Epileptic.or.Control', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1051 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0002629 0.00026288 2.0869 10000 0.1541
## Residuals 96 0.0120927 0.00012597
a table of # of breed. frequency table
(data.frame(sample_data(ps.prop))$Breed.Group..1.) %>% table()
## .
## Herder NA Pointer Spaniel Retriever Scent Hound
## 20 9 26 26 5
## Sight Hound Sled Dog Terrier
## 3 6 3
plot_ord(ps.prop, 'Breed.Group..1.','MDS','bray')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV546 -- in the phylogenetic tree in the data you provided.
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','NMDS','bray')
## Run 0 stress 0.2086595
## Run 1 stress 0.2098432
## Run 2 stress 0.2241091
## Run 3 stress 0.2274123
## Run 4 stress 0.2041237
## ... New best solution
## ... Procrustes: rmse 0.05286805 max resid 0.4152318
## Run 5 stress 0.2100706
## Run 6 stress 0.2039644
## ... New best solution
## ... Procrustes: rmse 0.03389528 max resid 0.2871853
## Run 7 stress 0.2092792
## Run 8 stress 0.4119602
## Run 9 stress 0.204043
## ... Procrustes: rmse 0.01527918 max resid 0.1203596
## Run 10 stress 0.2030238
## ... New best solution
## ... Procrustes: rmse 0.0153307 max resid 0.07956405
## Run 11 stress 0.2045171
## Run 12 stress 0.2247327
## Run 13 stress 0.2304353
## Run 14 stress 0.2041598
## Run 15 stress 0.2045553
## Run 16 stress 0.2139464
## Run 17 stress 0.2047252
## Run 18 stress 0.2102916
## Run 19 stress 0.2033599
## ... Procrustes: rmse 0.01554262 max resid 0.07812185
## Run 20 stress 0.2042144
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmax
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2224 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1709003
## Run 1 stress 0.1797137
## Run 2 stress 0.1709466
## ... Procrustes: rmse 0.005237464 max resid 0.03609124
## Run 3 stress 0.1740929
## Run 4 stress 0.1785395
## Run 5 stress 0.1960663
## Run 6 stress 0.1873395
## Run 7 stress 0.1755881
## Run 8 stress 0.1791391
## Run 9 stress 0.1711105
## ... Procrustes: rmse 0.01844282 max resid 0.1694661
## Run 10 stress 0.18166
## Run 11 stress 0.1871563
## Run 12 stress 0.174832
## Run 13 stress 0.1707072
## ... New best solution
## ... Procrustes: rmse 0.02223152 max resid 0.169232
## Run 14 stress 0.1837557
## Run 15 stress 0.1843321
## Run 16 stress 0.1788305
## Run 17 stress 0.1848787
## Run 18 stress 0.1748201
## Run 19 stress 0.1794211
## Run 20 stress 0.1794372
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmax
## Too few points to calculate an ellipse
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Too few points to calculate an ellipse
# same breed group v.s. diff breed group in household
# check if data are paired within household
if (!identical(sample_data(ps)$Household[seq(1,98,2)],
sample_data(ps)$Household[seq(1,98,2)+1])) stop('data is not paired')
same.breed.group <- table(sample_data(ps)$Household, sample_data(ps)$Breed.Group..1.) %>%
apply(1, function(x) any(x == 2)) # if dog's grouping are same with in each house
same.breed.group;sum(same.breed.group)
## 1 10 11 12 13 14 15 16 17 18 19 2 20
## FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
## 21 22 23 24 25 26 27 28 29 3 30 31 32
## TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE
## 33 34 35 36 38 39 4 40 41 42 43 44 45
## TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 47 48 49 5 50 51 6 7 8 9
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [1] 39
Here we see 39 out of 49 household have dogs in same breed group
permanova(ps.prop, 'Household + Breed.Group..1.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.3397 9.999e-05 ***
## Breed.Group..1. 6 0.8417 0.04746 1.2896 0.07649 .
## Residual 43 4.6775 0.26374
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV141 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 1.31342 0.67548 2.2590 9.999e-05 ***
## Breed.Group..1. 6 0.11016 0.05666 1.5158 0.047 *
## Residual 43 0.52084 0.26786
## Total 97 1.94442 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Here we
ps.control <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Control')
ps.epileptic <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Epileptic')
perm.breed.bray.ctl <- permanova(ps.control, 'Breed.Group..1.', 'bray')
perm.breed.wunif.ctl <- permanova(ps.control, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2095 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.ctl;perm.breed.wunif.ctl
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.7534 0.194 1.4098 0.0314 *
## Residual 41 7.2851 0.806
## Total 48 9.0385 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.16801 0.20757 1.5343 0.0409 *
## Residual 41 0.64140 0.79243
## Total 48 0.80942 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm.breed.bray.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'bray')
perm.breed.wunif.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV74 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.epi;perm.breed.wunif.epi
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.3363 0.15645 1.0863 0.2502
## Residual 41 7.2055 0.84355
## Total 48 8.5418 1.00000
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.56526 0.17885 1.2758 0.1108
## Residual 41 2.59518 0.82115
## Total 48 3.16044 1.00000
\[X^2=-2 \sum_{i=1}^k \ln \left(p_i\right)\]
bray.p <- c(perm.breed.bray.ctl$`Pr(>F)`[1], perm.breed.bray.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(bray.p)), df = 4, lower.tail = FALSE)
## [1] 0.04592367
wunif.p <- c(perm.breed.wunif.ctl$`Pr(>F)`[1], perm.breed.wunif.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(wunif.p)), df = 4, lower.tail = FALSE)
## [1] 0.02898295
With using the weighted-unifrac distance, there’s at
least one of the null hypothesis is rejected.
Here we are only focusing on epileptic dogs.
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','bray')
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1518 -- in the phylogenetic tree in the data you provided.
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','bray')
## Run 0 stress 0.2020935
## Run 1 stress 0.2040404
## Run 2 stress 0.2100444
## Run 3 stress 0.2003524
## ... New best solution
## ... Procrustes: rmse 0.02934314 max resid 0.1297959
## Run 4 stress 0.1999022
## ... New best solution
## ... Procrustes: rmse 0.06175294 max resid 0.3320622
## Run 5 stress 0.2112704
## Run 6 stress 0.2201524
## Run 7 stress 0.2186232
## Run 8 stress 0.2107734
## Run 9 stress 0.2205906
## Run 10 stress 0.2041888
## Run 11 stress 0.2082349
## Run 12 stress 0.1992016
## ... New best solution
## ... Procrustes: rmse 0.05964454 max resid 0.3435966
## Run 13 stress 0.2065769
## Run 14 stress 0.1994009
## ... Procrustes: rmse 0.06828334 max resid 0.3684565
## Run 15 stress 0.2000371
## Run 16 stress 0.2073667
## Run 17 stress 0.2051983
## Run 18 stress 0.2155364
## Run 19 stress 0.2041885
## Run 20 stress 0.2118414
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 2: no. of iterations >= maxit
## 18: stress ratio > sratmax
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV530 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1938893
## Run 1 stress 0.1886968
## ... New best solution
## ... Procrustes: rmse 0.07335568 max resid 0.2991985
## Run 2 stress 0.1940666
## Run 3 stress 0.1923904
## Run 4 stress 0.2041958
## Run 5 stress 0.2038335
## Run 6 stress 0.1838299
## ... New best solution
## ... Procrustes: rmse 0.07801956 max resid 0.3016581
## Run 7 stress 0.1915817
## Run 8 stress 0.194143
## Run 9 stress 0.1923897
## Run 10 stress 0.1902326
## Run 11 stress 0.1867653
## Run 12 stress 0.1840865
## ... Procrustes: rmse 0.06725676 max resid 0.3077652
## Run 13 stress 0.1942088
## Run 14 stress 0.1924326
## Run 15 stress 0.1866678
## Run 16 stress 0.1991867
## Run 17 stress 0.2027093
## Run 18 stress 0.1837112
## ... New best solution
## ... Procrustes: rmse 0.02928892 max resid 0.1736966
## Run 19 stress 0.1946872
## Run 20 stress 0.1843385
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.epileptic, 'Pheno.Y.N','bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.2211 0.02588 1.2488 0.1821
## Residual 47 8.3207 0.97412
## Total 48 8.5418 1.00000
permanova(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1824 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.000493 0.01497 0.7141 0.6772
## Residual 47 0.032422 0.98503
## Total 48 0.032914 1.00000
permdisp(ps.epileptic, 'Pheno.Y.N','bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.02029 0.020293 1.5629 10000 0.2212
## Residuals 47 0.61028 0.012985
permdisp(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2641 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01049 0.0104903 2.803 10000 0.1053
## Residuals 47 0.17590 0.0037426
# within epileptic dogs
table(sample_data(ps)$Seizure.Freedom..Y.N., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## No 0 32
## Yes 0 14
table(sample_data(ps)$Seizure.Control..Satisfactory.Unsatisfactory., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## S 0 33
## US 0 13
Does health condition related to dog’s sex? Here we conduct a Chi-squared test.
# chisq test on epi male dog v.s. female.
# not within health condition
tb <- table(sample_data(ps)$Sex, sample_data(ps)$Epileptic.or.Control)
tb
##
## Control Epileptic
## F 33 25
## M 16 24
chisq.test(tb)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tb
## X-squared = 2.0698, df = 1, p-value = 0.1502
plot_ord(ps, 'Sex','MDS','bray')
plot_ord(ps, 'Sex','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV762 -- in the phylogenetic tree in the data you provided.
plot_ord(ps, 'Sex','NMDS','bray')
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.242098
## Run 1 stress 0.2405582
## ... New best solution
## ... Procrustes: rmse 0.03689367 max resid 0.2477997
## Run 2 stress 0.2535572
## Run 3 stress 0.2506086
## Run 4 stress 0.2481361
## Run 5 stress 0.2435657
## Run 6 stress 0.2485509
## Run 7 stress 0.2605266
## Run 8 stress 0.250288
## Run 9 stress 0.2520829
## Run 10 stress 0.2475608
## Run 11 stress 0.2457253
## Run 12 stress 0.2547584
## Run 13 stress 0.2559922
## Run 14 stress 0.2443796
## Run 15 stress 0.2478656
## Run 16 stress 0.2503331
## Run 17 stress 0.247363
## Run 18 stress 0.2442408
## Run 19 stress 0.2468106
## Run 20 stress 0.2503928
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmax
plot_ord(ps, 'Sex','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV103 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1726779
## Run 1 stress 0.1878752
## Run 2 stress 0.191488
## Run 3 stress 0.1707686
## ... New best solution
## ... Procrustes: rmse 0.0243061 max resid 0.1685687
## Run 4 stress 0.1761118
## Run 5 stress 0.1804319
## Run 6 stress 0.1881514
## Run 7 stress 0.1799698
## Run 8 stress 0.1700946
## ... New best solution
## ... Procrustes: rmse 0.01749301 max resid 0.1696919
## Run 9 stress 0.1798839
## Run 10 stress 0.1700947
## ... Procrustes: rmse 9.881015e-05 max resid 0.0006521253
## ... Similar to previous best
## Run 11 stress 0.1884663
## Run 12 stress 0.4119596
## Run 13 stress 0.1942981
## Run 14 stress 0.1766625
## Run 15 stress 0.188327
## Run 16 stress 0.1855331
## Run 17 stress 0.178715
## Run 18 stress 0.1792163
## Run 19 stress 0.183729
## Run 20 stress 0.1738765
## *** Best solution repeated 1 times
permanova(ps.prop, 'Household + Sex', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.2373 9.999e-05 ***
## Sex 1 0.0589 0.00332 0.5174 0.9769
## Residual 48 5.4604 0.30787
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1940 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 0.057810 0.67238 2.0607 9.999e-05 ***
## Sex 1 0.000113 0.00132 0.1942 0.9939
## Residual 48 0.028054 0.32630
## Total 97 0.085977 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.prop, 'Sex', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01829 0.018289 1.2389 10000 0.2669
## Residuals 96 1.41725 0.014763
permdisp(ps.prop, 'Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1341 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00026 0.0002634 0.0607 10000 0.8127
## Residuals 96 0.41665 0.0043401
sam$Age..months. <- as.numeric(sam$Age..months.)
ggplot(sam) +
geom_line(aes(x = as.numeric(Household), y = Age..months./12, group = Household)) +
geom_point(aes(x = as.numeric(Household), y = Age..months./12, colour = Epileptic.or.Control)) +
xlab('Household') + ylab('Age in Year')
sam.epi <- sam %>% filter(Epileptic.or.Control == 'Epileptic')
sam.ctl <- sam %>% filter(Epileptic.or.Control == 'Control')
age.diff <- data.frame(Household = sam.epi$Household, Age.Diff = sam.epi$Age..months. - sam.ctl$Age..months.)
ggplot(age.diff) +
geom_bar(aes(x = as.numeric(Household), y = Age.Diff/12), stat = 'identity') +
xlab('Household') + ylab('Age Difference in Year')
permanova(ps.prop, 'Household + Age..months.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.3084 9.999e-05 ***
## Age..months. 19 2.2116 0.12470 1.0557 0.3173
## Residual 30 3.3076 0.18650
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Age..months.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1132 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 4.2025 0.68111 2.2374 9.999e-05 ***
## Age..months. 19 0.7936 0.12863 1.0675 0.3479
## Residual 30 1.1739 0.19026
## Total 97 6.1701 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1